sound file SchroederCallsTest4_3Nov2022.wav:
schroeders <- read_wave("./data/raw/SchroederCallsTest4_3Nov2022.wav")
spectro(schroeders, scale = FALSE, osc = TRUE, fastdisp = TRUE, grid = FALSE,
palette = viridis, collevels = seq(-100, 0, 5), wl = 512, colwave = wave_col,
heights = c(1, 1))Spectro + oscillogram
bout_1 <- cutw(schroeders, from = 0, to = 0.2, output = "Wave")
spectro(bout_1, scale = FALSE, osc = TRUE, fastdisp = TRUE, grid = FALSE,
palette = viridis, collevels = seq(-100, 0, 5), wl = 512/2, colwave = wave_col,
heights = c(1, 1))Oscillogram
Oscillogram zooming in
Try it on the first bout. First get correlation for an amplitude template with the first few samples
# make a template from the first fifth of duration of the entire
# shroeder
template <- bout_1@left[seq_len(length(bout_1@left)/5)]
cors <- vapply(seq_len(length(bout_1@left)/4), function(x) {
segment <- bout_1@left[x:(x + length(template) - 1)]
cor(template, segment)
}, FUN.VALUE = numeric(1))
plot(cors, type = "l", col = "orange", xlab = "samples", ylab = "Pearson correlation")Then get location of correlation peaks
# get highest peak
tpks <- tpks[tpks[, 2] > 0.8, ]
# order by sample number
tpks <- tpks[order(tpks[, 1]), ]
# remove first peak
tpks <- tpks[-1, ]Then get the mean difference (in samples) between peaks (which is inversely related to the frequency)
## [1] 244
Now get the start of each period in the sound clip
# get envelope peaks
amppk <- fpeaks(cbind(1:length(bout_1@left), abs(bout_1@left)/max(abs(bout_1@left))),
plot = FALSE)
# filter highest peaks
amppk <- amppk[amppk[, 2] > 0.95, ]
# get the highest one
first_max <- amppk[which.max(amppk[, 2]), 1]/bout_1@samp.rate
# plot
oscillo(bout_1, from = 0, to = 0.03, colwave = wave_col)
abline(v = first_max + (-100:100 * (mean(sample_diffs)/bout_1@samp.rate)),
col = "red", lty = 2, lwd = 2)spot_periods <- function(wave, plot = TRUE, abs.amp = TRUE, alpha = 0.5,
from = 0, to = 0.03) {
# make a template from the first fifth of duration of the
# entire shroeder
template <- wave@left[seq_len(length(wave@left)/5)]
cors <- vapply(seq_len(length(wave@left)/4), function(x) {
segment <- wave@left[x:(x + length(template) - 1)]
cor(template, segment)
}, FUN.VALUE = numeric(1))
tpks <- fpeaks(cbind(1:length(cors), cors), plot = FALSE)
# get highest peak
tpks <- tpks[tpks[, 2] > 0.8, ]
# order by sample number
tpks <- tpks[order(tpks[, 1]), ]
# remove first peak
tpks <- tpks[-1, ]
sample_diffs <- diff(tpks[, 1])
# mean distance in samples between amplitude peaks
mean_diff <- mean(sample_diffs)
norm_abs <- if (abs.amp)
abs(wave@left)/max(abs(wave@left)) else wave@left/max(wave@left)
# get envelope peaks
amppk <- fpeaks(cbind(1:length(wave@left), norm_abs), plot = FALSE)
# filter highest peaks
amppk <- amppk[amppk[, 2] > 0.95, ]
# get the highest one
first_max <- amppk[which.max(amppk[, 2]), 1]/wave@samp.rate
# get starts of periods
positions <- first_max + (-1000:1000 * (mean_diff/wave@samp.rate))
# filter positions within duration of wave
positions <- positions[positions > 0 & positions < duration(wave)]
# plot
if (plot) {
osc <- oscillo(wave, from = from, to = to, colwave = wave_col)
abline(v = positions, col = adjustcolor("blue", alpha.f = alpha),
lty = 2, lwd = 3)
points(y = wave@left, x = seq_len(length(wave))/wave@samp.rate,
type = "l", col = wave_col)
}
return(positions)
}Try function on other shroeders
5th bout
bout_5 <- cutw(schroeders, from = 4.8, to = 5, output = "Wave")
pos <- spot_periods(wave = bout_5, alpha = 0.5)9th bout
bout_9 <- cutw(schroeders, from = 10.2, to = 10.4, output = "Wave")
pos <- spot_periods(wave = bout_9, to = 0.02)11th bout
bout_11 <- cutw(schroeders, from = 15, to = 15.2, output = "Wave")
pos <- spot_periods(wave = bout_11, to = 0.02)20th bout
bout_20 <- cutw(schroeders, from = 31.8, to = 32, output = "Wave")
pos <- spot_periods(wave = bout_20, to = 0.015)# make a template from the first fifth of duration of the entire
# shroeder
template <- bout_1@left[seq_len(length(bout_1@left)/8)]
dists <- pbapply::pbsapply(cl = 4, seq_len(length(bout_1@left)/4),
function(x) {
segment <- bout_1@left[x:(x + length(template) - 1)]
dtw_dist <- dtw::dtwDist(mx = rbind(template, segment))
return(dtw_dist[1, 2])
})
saveRDS(dists, "./data/processed/dtw_distance_bout_1.RDS")dists <- readRDS("./data/processed/dtw_distance_bout_1.RDS")
dists <- dists/max(dists)
sims <- 1 - dists
plot(sims, type = "l", col = "orange", xlab = "samples", ylab = "DTW similarities")Then get location of correlation peaks
Making the time series stationary (diff())
stat_sims <- diff(sims)
stat_sims <- stat_sims + abs(min(stat_sims))
stat_sims <- stat_sims/max(stat_sims)
tpks <- fpeaks(cbind(1:length(stat_sims), stat_sims), threshold = 0.98)# get highest peak
tpks <- tpks[tpks[, 2] > 0.98, ]
sels <- data.frame(sound.files = 1, selec = 1:nrow(tpks), start = tpks[,
1], end = tpks[, 1] + 6, peak = tpks[, 2])
# merge peaks close to each other
sels <- ohun::merge_overlaps(sels, pb = FALSE)
sels$dur <- sels$end - sels$start
tpks <- cbind(sels$start + (sels$dur - 6)/2, sels$peak)
# order by sample number
tpks <- tpks[order(tpks[, 1]), ]
# remove first peak tpks <- tpks[-1, ]Then get the difference (in samples) between peaks (which is inversely related to the frequency)
Now get the start of each period in the sound clip
# # get envelope peaks amppk <-
# fpeaks(cbind(1:length(bout_1@left), abs(bout_1@left) /
# max(abs(bout_1@left))), plot = FALSE) # filter highest peaks
# amppk <- amppk[amppk[,2] > 0.95, ] # get the highest one
# first_max <- amppk[which.max(amppk[,2]),1] / bout_1@samp.rate
# plot
oscillo(bout_1, from = 0, to = 0.03, colwave = wave_col)
abline(v = cumsum(sample_diffs)/bout_1@samp.rate, col = "red", lty = 2,
lwd = 2)# plot
oscillo(bout_1, from = 0, to = 0.045, colwave = wave_col)
abline(v = cumsum(sample_diffs)/bout_1@samp.rate, col = "red", lty = 2,
lwd = 2)
Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_CR.UTF-8 LC_COLLATE=pt_BR.UTF-8
## [5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] viridis_0.6.2 viridisLite_0.4.1 warbleR_1.1.28
## [4] NatureSounds_1.0.4 seewave_2.2.0 tuneR_1.4.1
## [7] xaringanExtra_0.7.0 rprojroot_2.0.3 formatR_1.11
## [10] knitr_1.42 kableExtra_1.3.4
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.9 svglite_2.1.0 fftw_1.0-7 assertthat_0.2.1
## [5] digest_0.6.31 packrat_0.9.0 utf8_1.2.2 R6_2.5.1
## [9] Sim.DiffProc_4.8 signal_0.7-7 evaluate_0.20 highr_0.10
## [13] httr_1.4.4 ggplot2_3.4.0 pillar_1.8.1 rlang_1.0.6
## [17] uuid_1.1-0 rstudioapi_0.14 jquerylib_0.1.4 rmarkdown_2.20
## [21] webshot_0.5.4 sketchy_1.0.2 stringr_1.5.0 igraph_1.3.5
## [25] RCurl_1.98-1.9 munsell_0.5.0 proxy_0.4-27 Deriv_4.1.3
## [29] compiler_4.1.0 xfun_0.36 pkgconfig_2.0.3 systemfonts_1.0.4
## [33] htmltools_0.5.4 tidyselect_1.2.0 gridExtra_2.3 tibble_3.1.8
## [37] dtw_1.23-1 fansi_1.0.3 crayon_1.5.2 dplyr_1.0.10
## [41] MASS_7.3-54 bitops_1.0-7 grid_4.1.0 DBI_1.1.1
## [45] jsonlite_1.8.4 gtable_0.3.1 lifecycle_1.0.3 magrittr_2.0.3
## [49] scales_1.2.1 cli_3.6.0 stringi_1.7.12 cachem_1.0.6
## [53] pbapply_1.6-0 remotes_2.4.2 xml2_1.3.3 bslib_0.4.2
## [57] vctrs_0.5.2 generics_0.1.3 ohun_0.1.0 rjson_0.2.21
## [61] tools_4.1.0 glue_1.6.2 parallel_4.1.0 fastmap_1.1.0
## [65] yaml_2.3.7 colorspace_2.0-3 rvest_1.0.3 sass_0.4.5